#ifndef CPFUNCTIONS_H
#define CPFUNCTIONS_H

#include <TObject.h>
#include <TMath.h>
class TF3;

class CpFunctions : public TObject {
public:
   CpFunctions();
   ~CpFunctions();

   static const enum EFunctionType { kReal=0,kRealDimLess=1, kTest=2} x;

	static void Init();
   static void PrintMagneticInfo(Double_t x,Double_t y,Double_t z);

	// constans
	static Double_t kMuZero;
	static Double_t kIs;
	static Double_t kParticleRho;
	static Double_t kParticleR;
	static Double_t kBloodViscosity;
	static Double_t kBloodVelX;
	static Double_t kBloodVelY;
	static Double_t kBloodVelZ;

  static Double_t kBulharConstant;
  static Double_t kMultipier;
  static Bool_t fgIsInMaget;

  

	static Double_t ParticleMass() { return 4*TMath::Pi()/3*TMath::Power(kParticleR,3)*kParticleRho;}
	static Double_t Chi(Double_t x,Double_t y,Double_t z);
	static TF3 *fgMagForceX;
	static TF3 *fgMagForceY;
	static TF3 *fgMagForceZ;

	// 1st derivates
	static TF3 *fgMagForceXdX;
	static TF3 *fgMagForceXdY;
	static TF3 *fgMagForceXdZ;
	static TF3 *fgMagForceYdX;
	static TF3 *fgMagForceYdY;
	static TF3 *fgMagForceYdZ;
	static TF3 *fgMagForceZdX;
	static TF3 *fgMagForceZdY;
	static TF3 *fgMagForceZdZ;

	// 2nd derivates (only needed)
	static TF3 *fgMagForceXdX2;
	static TF3 *fgMagForceXdY2;
	static TF3 *fgMagForceXdZ2;
	static TF3 *fgMagForceYdX2;
	static TF3 *fgMagForceYdY2;
	static TF3 *fgMagForceYdZ2;
	static TF3 *fgMagForceZdX2;
	static TF3 *fgMagForceZdY2;
	static TF3 *fgMagForceZdZ2;


	
	
	static Double_t dIp0(Double_t t,Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);
	static Double_t dIp1(Double_t t,Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);
	static Double_t dIp2(Double_t t,Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);
	static Double_t dIp3(Double_t t,Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);
	static Double_t dIp4(Double_t t,Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);
	static Double_t dIp5(Double_t t,Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);


	static Double_t dF0D0(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF0D1(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF0D2(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF0D3(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 1.0; }
	static Double_t dF0D4(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF0D5(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }


	static Double_t dF1D0(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF1D1(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF1D2(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF1D3(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF1D4(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 1.0; }
	static Double_t dF1D5(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }

	static Double_t dF2D0(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF2D1(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF2D2(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF2D3(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF2D4(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF2D5(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 1.0; }


	static Double_t dF3D0(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF3D1(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF3D2(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return -6.0*TMath::Pi()*kMuZero*kParticleR/ParticleMass();} // -6.0*Pi*nu*Radius/mp
	static Double_t dF3D3(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF3D4(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF3D5(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }

	static Double_t dF4D0(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF4D1(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF4D2(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF4D3(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF4D4(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return -6.0*TMath::Pi()*kMuZero*kParticleR/ParticleMass();} // -6.0*Pi*nu*Radius/mp
	static Double_t dF4D5(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }

	static Double_t dF5D0(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF5D1(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF5D2(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF5D3(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF5D4(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t dF5D5(Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return -6.0*TMath::Pi()*kMuZero*kParticleR/ParticleMass();} // -6.0*Pi*nu*Radius/mp


	//  function 6:   df_i/dt
	static Double_t df0dt(Double_t /*t*/, Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t df1dt(Double_t /*t*/, Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t df2dt(Double_t /*t*/, Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t df3dt(Double_t /*t*/, Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t df4dt(Double_t /*t*/, Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }
	static Double_t df5dt(Double_t /*t*/, Double_t /*x*/,Double_t /*y*/,Double_t /*z*/,Double_t /*vx*/,Double_t /*vy*/,Double_t /*vz*/) { return 0.0; }

   static Double_t MyNabla1(Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);
   static Double_t MyNabla2(Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);
   static Double_t MyNabla3(Double_t x,Double_t y,Double_t z,Double_t vx,Double_t vy,Double_t vz);

   static Double_t BloodVelocityAverage();

	static int GslFunc(double t, const double Ip[], double f[], void *params);
	static int GslJac(double t, const double Ip[], double *dfdy, double dfdt[], void *params);

	static void SetType(EFunctionType t) {fgType = t;}

private:
	static EFunctionType fgType;

//    ClassDef(CpFunctions, 1);
};

#endif
